function h = mcgh_histfit(data,nbins,mr1,sr1,p1,mr0,sr0,p0,mr_1,sr_1,p_1)
%HISTFIT Histogram with superimposed fitted normal density.
%   HISTFIT(DATA,NBINS) plots a histogram of the values in the vector DATA.
%   using NBINS bars in the histogram. With one input argument, NBINS is set 
%   to the square root of the number of elements in DATA. 
%
%   H = HISTFIT(DATA,NBINS) returns a vector of handles to the plotted lines.
%   H(1) is a handle to the histogram, H(2) is a handle to the density curve.
%   B.A. Jones 2-14-95
%   Copyright 1993-2000 The MathWorks, Inc. 
%   $Revision: 2.12 $  $Date: 2000/05/26 18:52:57 $
global p_1 p1 p0 m1 m0 m_1 s_1 s1 s0
global numofrun

[row,col] = size(data);
if min(row,col) > 1
   error('First argument has to be a vector.');
end

if row == 1
  data = data(:);
end
row = sum(~isnan(data));

if nargin < 2
  nbins = ceil(sqrt(row));
end

%add wang 28/01/2002
[n,xbin]=mcgh_hist(data,nbins);
%end add wang

mr = mcgh_nanmean2(data); % Estimates the parameter, MU, of the normal distribution.
sr = mcgh_nanstd2(data);  % Estimates the parameter, SIGMA, of the normal distribution.

%x=(-3*sr+mr:0.1*sr:3*sr+mr)';% Evenly spaced samples of the expected data range 3.
%default is above function with 3std from mean , I changed to new range
x=(min(min(data)):0.1*sr:max(max(data)))';
intp_n=interp1(xbin,n,x,'nearest');
intp_n(find(isnan(intp_n)==1))=0; %replace nan as zero

 hh = bar(xbin,n,1); % Plots the histogram.  No gap between bars.
 np = get(gca,'NextPlot');    
 set(gca,'NextPlot','add')    
                             
 xd = get(hh,'Xdata'); % Gets the x-data of the bins.

 rangex = max(xd(:)) - min(xd(:)); % Finds the range of this data.
 binwidth = rangex/nbins;    % Finds the width of each bin.

%add aug 17
%disp('Start EM.........')
%numofrun
if numofrun<1
    [xx,yy]=mcgh_EMforGaussian(data,binwidth,row);
    % hh1 = plot(xx,yy,'r-','LineWidth',2);     % Plots density line over histogram.
    hold on
    plot(x,intp_n,'g-.','LineWidth',1.5);
    
    idx_yy=find(yy<=max(intp_n));
    hold on;
    hh1 = plot(xx(idx_yy),yy(idx_yy),'r-','LineWidth',2);     % Plots density line over histogram.
    disp('EM estimate');
else
   

    %added wang 28/01/2002
    %mr1=4;
    %sr1=0.2;
    %   p1=0.2
    %mr0=0;
    %sr0=0.2;
    %p0=0.8
    %mr_1=-1;
    %sr_1=0.2;
    %p_1=0.0
    %   y = normpdf(x,mr,sr);  
    y=p1.*mcgh_normpdf2(x,mr1,sr1)+p0.*mcgh_normpdf2(x,mr0,sr0)+p_1.*mcgh_normpdf2(x,mr_1,sr_1);  
    y = row*(y*binwidth);   % Normalization necessary to overplot the histogram.
    %end added

    %added wang 28/01/2002
    y=y./row;
    %end added

    %hh1 = plot(x,y,'r-','LineWidth',2);     % Plots density line over histogram.
    %[x,y]
    idx_yy=find(y<=max(intp_n));
    hold on;
    hh1 = plot(x(idx_yy),y(idx_yy),'r-','LineWidth',2);     % Plots density line over histogram.
    %disp('EM estimate');
    
    hold on
    plot(x,intp_n,'g-.','LineWidth',1.5);
    %[x,intp_n]
end %end if

if nargout == 1
  h = [hh; hh1];
end

set(gca,'NextPlot',np) 
xl=xlabel('Log2 ratio');
yl=ylabel('Ratio distribution');
set(xl,'FontSize',13);
set(yl,'FontSize',13);
